#!/usr/bin/perl -w
use strict;

# run like this where the prefix of the output file will be tter_
# and the RepeatMasker out file is Thielavia_terrestris_NRRL-8126_1.fasta.out
# and the FASTA genome file is Thielavia_terrestris_NRRL-8126_1.fasta

# perl extract_RepeatMasker_features.pl -db Thielavia_terrestris_NRRL-8126_1.fasta -p tter_ Thielavia_terrestris_NRRL-8126_1.fasta.out

use Bio::DB::Fasta;
use Bio::SeqIO;
use Getopt::Long;

my $db;
my $prefix = './';
GetOptions(
	   'd|db:s' => \$db,
	   'p|prefix:s' => \$prefix,
	   );

if( ! defined $db ) {
    die"Must provide a genome db with -d or --db option\n";
}
my $dbh = Bio::DB::Fasta->new($db);
my %fams;
my $n =0;
while(<>) {
    next if /^\s*(SW|score)/ || /^\s+$/;
    my ($score,$div,$del,$ins,$query,$qstart,$qend,$q_left,$strand,
	$family, $class, $hstart,$hend,$h_left,$id) = split;
    push @{$fams{$class}},[$family,$query,$qstart,$qend,$strand];
    $n++;
}

for my $class ( keys %fams ) {
    print "$class ",scalar @{$fams{$class}},"\n";
    my $classid = $class;
    $classid =~ s/\//_/; # let's change the / in LINE/Tad1 to LINE_Tad1
    my $out = Bio::SeqIO->new(-format => 'fasta',
			      -file   => ">$prefix$classid.fas");
    for my $location ( @{$fams{$class}} ) {
	my ($family,$chrom,$start,$end,$strand) = @$location;
	if( $strand eq 'C' ) {
	    # flip start/end if on opposite strand	
	    ($start,$end) = ($end,$start);
	}
	
	# get the sequence from the genome db
	my $seq = $dbh->seq($chrom, $start => $end);
	my $seqobject = Bio::PrimarySeq->new(-id => "$chrom\_$start\_$end",
					     -desc => $family,
					     -seq => $seq);
	$out->write_seq($seqobject);
    }
}
